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ABSTRACT 

Context. With the plethora of detailed results from heliospheric missions such as Ulysses and SOHO and at the advent of the first 
mission dedicated to in situ studies of neutral heliospheric atoms IBEX, we have entered the era of precision heliospheric studies. 
Interpretation of these data require precision modeling, with second-order effects quantitatively taken into account. 
Aims. We study the influence of the non-flat shape of the solar Lyman-a line on the distribution of neutral interstellar hydrogen in the 
inner heliosphere and assess the importance of this effect for interpretation of heliospheric in situ measurements. 
Methods. Based on available data, we (i) constructed a model of evolution for the solar Lyman-fy line profile with solar activity, (ii) 
modified an existing test-particle code used to calculate the distribution of neutral interstellar hydrogen in the inner heliosphere so that 
it takes the dependence of radiation pressure on radial velocity into account, and (iii) compared results of the old and new version. 
Results. Discrepancies between the classical and Doppler models appear between ~ 5 and 3 AU and increase towards the Sun from a 
few percent to a factor of 1.5 at 1 AU. The classical model overestimates the density everywhere except for a ~ 60° cone around the 
downwind direction, where a density deficit appears. The magnitude of the discrepancies appreciably depends on the phase of the solar 
cycle, but only weakly on the parameters of the gas at the termination shock. For in situ measurements of neutral atoms performed at 
~ 1 AU, like those planned for IBEX, the Doppler correction will need to be taken into account, because the modifications include both 
the magnitude and direction of the local flux by a few km/s and degrees, respectively, which, when unaccounted for, would introduce 
an error of a few km/s and degrees in determination of the magnitude and direction of the bulk velocity vector at the termination 
shock. 

Conclusions. The Doppler correction is appreciable for in situ observations of neutral H populations and their derivatives performed 
a few AU from the Sun. 

Key words. Sun: UV radiation - Ultraviolet: solar system -Interplanetary medium - Line: profiles - ISM: atoms 



1. Introduction 

Development of models of the distribution of neutral inter- 
stellar hydrogen in the heliosphere began after the discovery 
of a diffuse interplanetary L yman-o' glow (Thomas & Krassal 
197lHBert aux & Bl amontjl97lh . interpreted bv lFahjd 1968b and 
BlumjfcFahjI([i970) as caused by the scattering of solar Lyman- 
a radiation on this gas flowing through the Solar System. 

At the very early phase, the influence of processes going 
on at the heliospheric interface on the heliospheric gas was ne- 
glected. It was assumed that both the ionization due to solar out- 
put and the radiation pressure acting on the H atoms are sta- 
tionary and spherically symmetric around the Sun and that their 
rates fall off proportionally to inverse square of heliocentric dis- 
tance. It was further assumed that the inflowing neutral gas is 
monoenergetic, i.e., that before the encounter with the Sun the 
atoms move with identical, parallel-oriented velocities, equal to 
the macroscopic bulk velocity of the gas. 

These assumptions formed the basis of the first model of 
density distribution of neutral inte rstellar gas around the S un: 
the purely analytical "cold" model dFahril 1 9681: 1 Axfordll 1 9721) . It 
features axial symmetry about the axis of inflow of the gas and 
shows either a singularity at the downwind axis, when radiation 
pressure is too weak to compensate for solar gravity, or an empty 
"avoidance zone" in the downwind region with a paraboloidal 
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boundary surface, when radiation pressure overcompensates for 
solar gravity. 

Lifting of the monoenergetic assumption (i.e. allowing for a 
finite temperature of interstellar gas, high enough to yield ther- 
mal velocity comparable to the bulk velocity of the gas), adop- 
tion of the distribution function of the gas far away from the Sun 
("in infinity") in the form of a Maxwellian shifted in the velocity 
space by the bulk velocity vector, and assumption that the gas 
is collisionless on the distance scale comparable to the size of 
the heliosphere, all allowed us to use the Boltzmann equation to 
describe the problem of the interaction of neutral interstellar gas 
with the solar environment. Its solution brought t he "hot model" 
of the gas distribution in the inner he l iosphere dThoma s| 1978b 
lFahn[T978l fl979t IWu & Juded [19791: lLallement et alJll985al) . 
This model required numerical integration of the distribution 
function, but the function itself was still analytical. The "hot 
model" in its various implementations became the canonical 
model of neutral interstellar gas distribution in the inner he- 
liosphere, so further on we will refer to it as "the classical hot 
model". 

Further modeling of the interaction and distribution of neu- 
tral interstellar hydrogen near the Sun focused on two main top- 
ics. On one hand, a lot of effort was made toward simulating 
and understanding processes going on at the boundary between 
the expanding solar wind and the incoming, partially ionized in- 
terstellar gas, in the region referred to as the heliospheric inter- 
face - and this issue will not be addressed in the present paper: 
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the reader is refer r ed to recent reviews by Baranovl d2006blal) : 
llzmodenovl d2006l) ; llzmodenov & Baranovl (120061) . On the other 
hand, development of the hot model continued, aimed at a more 
realistic description of the distribution of neutral interstellar hy- 
drogen in the inner heliosphere, suitable for quantitative, and not 
only qualitative interpretation of heliospheric measurements. 

In the first shot, the assumption that the i onization rate is 
spheri cally symmteric was eliminated, when lLallement et al.l 
(1985b) described the latitudinal modulation of the charge ex- 
change rate with a one-parameter formula 1 - A sin 2 0, imple- 
menting it in the CNRS model of the heliospheric gas distribu- 
tion. This allowed the authors to vary the equator-to-pole con- 
trast of the ionization rate, but required keeping the width and 
range of the enhanced ionization band fixed. 

A different extension of t he hot model was proposed by 
IRucinski & Fahrl dl989l Il99ll) . who realized that the rate of 
ionization by electron impact is not proportional to the in- 
verse square of solar distance. Electro n ionization is particu- 
larly important for interstellar he lium dMcMullin et alJ 120041: 
Lalle ment et alJl2004t IWittel l2004'). The effects in hydro gen are 
only noticeable within a few AU from the Sun, where the den- 
sity of hydrogen gas is already very reduced by earlier ionization 
and solar radiation pressure, overcompensating for solar gravity; 
therefore this aspect of the physics of neutral interstellar hydro- 
gen has been neglected for quite a while since th en, but is rein- 
trodu ced in the most recent versions of the model (Bzowski et al. 
120081) . 

In the next round in the development of the density model, 
gone was assumption of the invariability of radiation pres- 
sure and ionization rate. This aspect of the physics of neu- 
tral interstellar gas was im p leme nted in a simplified numer- 
ical model by iKvrola et alJ d 19941) and shortly afterwards in 
the Warsaw test-particle code, to which we will refer in 
this paper as the Warsaw m od el fRuciriski & Bzowski [l995a; 
IBzowski & Rucinskil Il995allbl IBzowski el : al.l Il997l) . Indeed, 
both the solar EUV flux dFlovd et al.l l2002) and the solar wind 
flux dKing & P apitashvili 2005J) vary considerably during the so- 
lar cycle. The result is a solar cycle variation of the solar radi- 
ation pressure, of the EUV ionization rate, and of the rate of 
charge exchange between neutral H atoms and solar wind pro- 
tons. This in turn results in appreciable variations in the density 
and bulk velocity of neutral interstellar hydrogen within a dozen 
AU from the Sun. 

The Warsaw time-dependent model was fully numerical: not 
only did the integration of the distribution function have to be 
performed numerically, as in the classical "hot model", but also 
the calculation of the integrand, i.e. of the distribution function 
itself, was no longer analytical. At this phase in the heliospheric 
research, substitute idealized models of the evolution of these 
parameters had to be used due to the lack of sufficiently long 
time series of measurements of the solar wind speed and density, 
and of the solar EUV output. These aspects of the development 
of heliosp heric gas models were discu ssed in greater detail in a 
review by Rucihski & Bzowsk| d!996l) . 

In the next move, the inner-heliospheric effects of the inter- 
action of the solar wind and interstellar gas within the helio- 
spheric interface were taken into account. Owing to the charge 
exchange between the atoms of interstellar gas and the heated 
and compressed plasma in front of the heliopause, the origi- 
nal population of neutral atoms is somewhat cooled and ac- 
celerated, and a new population of atoms appears. These new 
neutral atoms initially maintain the properties of the parent 
plasma population, but later on they decouple from this plasma 
and flow through the heliopause and inside the termination 



shock of the solar wind. The two neutral populations interact 
further with the ionized components, exchanging charge with 
the protons from the local plasma and producing new pop- 
ulations of neutral atoms. Hence, the processes in the helio- 
spheric interface cr eate a few distinct, col l isionless populations 
of neutral atoms dOsterbart & Fahrlll992t iBaranov et ail 1 1991 ' 



Baran ov~& Malamall 1993b, as di scussed in detail by llzmodenovl 
d2000h and lMalama et alJ d2006l) . 

From the viewpoint of modeling the distribution of neu- 
tral interstellar hydrogen in the inner heliosphere, the most im- 
portant aspect of the processes going on in the heliospheric 
interface is the modifi cation of the di stribution function at 
the termination shock (Izmodeno v et alj|2001l) . Instead of the 
shifted Maxwellian with parameters homogeneous in space, 
IScherer et al.l d 19991) adopted a sum of two Maxwellians, with 
non-isotropic temperatures, shifted by appropriate bulk velocity 
vectors. The two components of the new functions corresponded 
to the two thermal populations (primary and secondary) pre- 
dicted by the Moscow Monte Carlo simulation of heliospheric 
interface and had parameters that depended on the offset angle 
9 from the upwind direction. This new version of the Warsaw 
model was time dependent, but still axially symmetric. At that 
time, enough measurements had been published to attempt to in- 
troduce observation-based models of the radiation pressure and 
ionization rate to the simulations. 

Axial symmetry was removed in the next round of the 
Warsaw model developm ent, when an an i sotropy of t he ioniza- 
tion rate was introduced. Bzows ki et al. | d200lL 12002) allowed 
the ionization rate to change as a continuous function of the 
heliographic latitude, with the latitudinal profile of the ioniza- 
tion rate continuously changing with the phase of solar cycle. 
As a result, the ionization field in the model became 2D (main- 
tainining the axial symmetry about the solar rotation axis) and 
time-dependent; and since the gas inflow axis and solar rotation 
axis are inclined at an angle to each other, the model of neutral 
hydrogen distribution in the inner heliosphere became 3D and 
time-dependent. 

The most recent extension of the Warsaw model is presented 
in this paper. We improve on the modeling of radiation pressure 
acting on individual H atoms, which now is not only a function 
of time, but also of the radial velocity of the atom. In the previous 
versions of the model it was assumed that the profile of the solar 
Lyman-ff line, responsible for the radiation pressure, is flat. In 
reality it not only is not flat, but shows co nsiderable variations 
depending on the p hase of solar activity dLemaire et alj [2002; 
iLemaire et aT]|2005l) . 

We take this into account, and in Sect. 2 we develop an 
observation-based model of the evolution of the line profile as 
a function of the wavelength- and disk-integrated flux. We use 
this model in a newly-developed code that simulates the density 
and higher moments of the distribution function of neutral inter- 
stellar hydrogen in the inner heliosphere. We discuss these cal- 
culations in Sect. 3. With the use of the newly-developed code, 
we assess modifications of the local gas density and flux with 
respect to the results of models neglecting the dependence of ra- 
diation pressure on the radial velocity. We also show possible 
implications for interpretation of heliospheric in situ measure- 
ments of neutral atoms and pickup ions. This discussion is pro- 
vided in Sect. 4. Section 5 offers a summary of the results. 



Tarnopolski & Bzowski: Effects of wavelength-dependent radiation pressure on heliospheric hydrogen 



3 



H solar Ly-a profiles 



of time, 



1.5 



0.5 



-200 -100 100 
radial velocity [km/s] 



200 



Fig. 1. Selected profiles of the solar Lyman-o' line observed by 
SUMER/SOHO (dots) for solar minimum (lower line, Feb. 17, 
1997) and solar maximum (upper line, May 20, 2000), compared 
with the results of the fitted model, specified in Eq. ((5]) (lines). 
The gray part of the model lines corresponds to +30 km/s around 
and illustrates the range on the profile that is the most rele- 
vant for the thermal populations of neutral H atoms in the he- 
liosphere. The horizontal lines mark the ft values used in the 
comparison simulations with the use of the classical hot model. 

H solar Ly-a profiles 



1.5 
1 

0.5 








!/% 








4 \ 

1 V 








1 1 

1 

1 













-200 -100 100 
radial velocity [km/s] 



200 



Fig. 2. Comparison of models of the solar Lyman-o' line profile: 
solid line: the present model; dashed line: Chabrillat & Kockart; 
dashed-dot line: Fahr et al.; dots: Scherer. 



2. The Doppler model of neutral hydrogen 
distribution in the inner heliosphere 

The approach exercised in t he present model is a modific ation 
of the approach presented by Ruci riski & BzowskH ([1995 a). The 
density and higher moments of the distribution function of neu- 
tral interstellar hydrogen at a location R in the inner heliosphere 
at a time r are calculated by numerical integration of the local 
distribution function, constructed as a product of the distribution 
function in the source region w stc and of the probability w; on of 
survival of the atoms traveling from the source region to the local 
point R against ionization, 

/ (R, V, t) = w src (r src (R, V, t) , v src (R, V, r)) w ion (R, V, r) , 

(1) 

where r SIC (R, V, t) is the start position in the source region of 
the atom that at the local point R at time t has velocity V, and 
V SIC (R, V, t) is the relevant start velocity in the source region. 

From the numerical point of view, the distribution function 
in the source region w STC can a priori be any reasonable function. 
In the following we adopted it either as a Maxwellian shifted by 
the bulk velocity vector, homogeneous in space and independent 
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or, in the two-populations case, as a sum of two such func- 
tions with different parameters n SIC (density in the source re- 
gion), vb (bulk velocity at the source region), and r src (temper- 
ature at the source region); also, m is the atom mass, and k the 
Boltzmann constant. Basically, the parameters of these popula- 
tions are functions of t he location in the heliospheric interface 
dlzmodenov et al. 2001). Since our goal in this paper is to inves- 
tigate the effect of the v r -dependence of radiation pressure on the 
gas distribution, we left this problem out of the analysis. Some 
insight into modifications of the gas density profile due to the 
angular gradient of the paramet ers of the H populatio ns at the 
termination shock is provided bv lBzowski et al.l d2008). 

The challenging issue in calculating the local distribution 
function of the gas is finding the link between the local veloc- 
ity of the test atom V at time r and location R and its location 
r src and velocity v slc in the source region in a scenario where 
the radiation pressure is a function of time and radial velocity. 
Indeed, the force acting on an H atom in the heliosphere is com- 
posed of the attracting solar gravity and the repelling radiation 
pressure, which is a function of the spectral flux lx relevant for a 
given radial velocity v r . The spectral flux lx is also a function of 
the line-integrated flux 7 tot (t), which varies on time scales from 
days to dozens of years. Thus the spectral flux responsible for the 
instantaneous radiation pressure acting on an H atom, moving 
with radial velocity v r , is a function of v r , / tot (f) and, implicitly, 
of t: I A = h (v r , / to t (*))• 

Since both / tot and solar gravity fall off with the heliocentric 
distance as 1/r 2 , the radiation pressure can be expressed by a fac- 
tor n of compensation of the solar gravity. In the modeling, this 
factor is a function of the radial velocity and line-integrated flux, 
which in turn is a function of time; hence, finally yi — p (y r , t). 

The link between the local position and velocity vectors R, V 
and the corresponding position and velocity vectors at the source 
region r src , v slc is determined by solving the equation of motion 
specified in Eq.©, performed numerically backwards in time. 
Along with the calculation of the trajectory of the atom, its prob- 
ability of survival is calculated. 
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In the equation above, J3o(r, t, <p) = r 2 f)(r, t, 0) is the scaled ion- 
ization rate at a heliographic latitude </>, r 2 = x 2 + y 2 + z 2 , 
GM is the constant of gravity times the solar mass, and v r (f) = 
(r(t) ■ v(t))/\r(t)\ is the radial velocity of the atom. The position 
vector r(f) has coordinates x, y, z and the velocity vector v(f) co- 
ordinates v x , v y , v z ; up is used to calculate the survival probability 
as: 

w ion = exp(-M /8 ). (4) 

Once the position and velocity vectors of the test atom in the 
source region are found (respectively, r src , v src ), the value of the 
distribution function in the source region w src can be calculated. 
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The value of the local distribution function is then obtained as 
a product of w src and the survival probability Wi on . To obtain the 
density and higher moments, the local distribution function is 
numerally integrated in the velocity space d 3 V. 

The model of radiation pressure (Tarnopolski 2007) was de- 
veloped based on a series of observations of the solar Lyman- 
a line profiles, performed by SUMER/SOHO between solar 
minimum and maximum dLemaire et alj [2002) and available 
on the Web. Since the SUMER instrument is located at the 
SOHO spacecraft orbiting aroun d the LI Lagrange point be- 
tween the Earth and the Sun dWilhelm et all 1 1 19991) . the ob- 
servations are free of the geocoronal contamination that af- 
fected earlier measurements carri ed out from low-Earth orbits 
(e.g. OSO-5, IVidal-Madia ri| l975l) The ab solute calibration of 
the wavelength reported by Lemaire et al.l (2002) is better than 
+0.0015 nm. The calibration of intensity was performed by di- 
rect comparison of results of integrating the observed profiles 
with the absolute fluxes obtained with the SOLSTICE instru- 
ment. The accuracy of the net flux was reported at a +10% 
level. It seems, however, that the data published on the Web 
are a smoothed version of actual measurements because the data 
points within the profiles do not show any scatter in spite of the 
accuracy of 10%, reported in the paper. 

Fittin g the func t ional forms of the profiles used in earlier 
studie s dFahj|1978b IChabrillat & Kockartsl [19971 IScherer et alJ 
2000) did not yield sufficient accuracy of the model (cf Fig.[T]and 
[2]). To develop the model of fi(v r , I lot ) discussed in this paper, the 
original data were rescaled to the spectral flux in the fi units as a 
function of radial velocity in km/s. A satisfactory re sult could be 
obtain ed with a 9-parameter function in the form (Tarnopolski 
l2007h 

/i (v,, /«(*)) = A [1 + B /„(»)] exp (-Cv 2 ) (5) 

x [l + D exp [Fv r - Gvj) + H exp i-Pv r - Qv 2 r )j 

with the following parameters: 

A = 2.4543 x 10~ 9 , B = 4.5694 x 10~ 4 , C = 3.8312 x 10~ 5 , 
D = 0.73879, F = 4.0396 x 10~ 2 , G = 3.5135 x 10~ 4 , 
# = 0.47817, P = 4.6841 x 10~ 2 , Q = 3.3373 x 10~ 4 . 
An interesting aspect of the model is that its sole solar- 
related parameter is the line-integrated flux in the Lyman-Q' line 
/tot- Any dependence of the radiation pressure on time goes to 
the model via / tot . Hence, with a model o f the behavior o f / to t 
in time in hand (such as, e.g., discussed bv lBzowskH d2001l) and 
iBzowski et al.l d2008l) ) we can immediately construct a model of 
the radiation pressure that will depend on the radial velocity of 
the atoms as well. 

The function defined in Eq.© with the parameters as listed 
above reproduces the observed profiles well both for solar min- 
imum (the profiles from July 27, 1996 to August 24, 1997) 
and solar maximum (the profiles from August 20, 1999 to 
August 22, 2001). The accuracy of the fit for the most inter- 
esting region ±140 km/s about the line center exceeds the ac- 
curacy of observations declared by the authors. Furthermore, 
the line-integrated fluxes obtained from the model profiles agree 
wit h the line-integrated data and with the total fluxes reported 
by iTobiska et al.l d2000t) for a given day. The differences are 
around 5%. An exception is the profile observed on May 20, 
2000, for which the difference is 9%, still within the abso- 
lute calibration accuracy of 10%. A comparison of the experi- 
mental data with the fitted model is shown i n Fig. Q] A com- 
parison of this model with earlier models (|Fahr et al.l I198U 
IChabrillat & Kockartslll997l:IScherer et alj|2000l) is presented in 

Rg. m 



Since the absolute calibration of / tot has changed since the 
time of publication of the profiles, we show the original profiles 
recalibrated by multiplying by an appropriate factor so that their 
absolute fluxes agree with the values accepted nowadays. The 
differences between the models in the velocity range +100 km/s 
around the line center are around 10 to 15%. For higher veloc- 
ities, the percentage differences between the models are not as 
important because the absolute values of radiation pressure are 
low anyway; furthermore, the fast atoms are not v ery sensitive 
to the radiation pressure in general (Bzowski 2008). The highest 
differences in the absolute terms occur at the peaks of the pro- 
file, i.e., at the velocities characteristic of the fast wing of the 
distribution function of thermal populations of neutral hydrogen 
in the inner heliosphere. With as many as 9 data sets available, 
covering the entire span of solar activity levels, we were able to 
come up with a model that seems to reproduce physical reality 
better than the former ones. 

Our model does not feature any nonlinear response of the 
line profile to / to t, although such a char acter of the response 
migh t be expected based on the findings of Ambl ard et al.l (f.2008) 
and lEmerich et alj (12005). Amblard et al. show that the spectrum 
of the Sun in the EUV region is composed of three components 
whose contributions to the net sprectral flux vary with the phase 
of the solar cycle, and they try to attribute these components 
to different features at the solar surface and in the solar atmo- 
sphere. Since the distribution of these features is a function of 
heliolatitude, one might then expect that the shape of the solar 
Lyman-a profile should vary with heliolatitude, but when lack- 
ing actual data, a construction of such a model at present is im- 
possible. Emerich et al. presented a nonlinear fit of the central 
value of the solar spectral flux in the Lyman-a line as a func- 
tion of the net flux / to t, but within the / tot range of 3.3 x 10 u to 
6.5 X 10 11 crrT 2 s _1 , deviations in their formula from a purely- 
linear dependence do not exceed 2%, which is much less than all 
the uncertainties of the profile and / tot measurements. Hence we 
believe the model we use in this paper is a reasonable approxi- 
mation, given the accuracy of available observations. 

The numerical scheme presented in this section was devised 
to enable calculating the density and higher moments of the dis- 
tribution function in the case of a fully 3D and time-dependent 
radiation pressure and ion ization rate. The code ha s alread y been 
used in its full mode by Tarnop olski & Bzowskil d2008l) to as- 
sess the flux of inter stellar deuterium at the Earth orbit and by 
Bzowski et al. (2008) to investigate the density of neutral inter- 
stellar H at the termination shock based on the observations of 
hydrogen pickup ions by Ulysses. In the present study, the code 
was restricted to the spherically symmetric and stationary case to 
facilitate an assessment of the influence of the v,.-dependence of 
the radiation pressure on the distribution of interstellar hydrogen 
in the inner heliosphere, unblended with other departures from 
the classical hot model. 



3. Calculations 

The simulations were performed in two groups. In the first one, 
we investigated the differences between the classical hot model 
and the Doppler model in a scenario with two populations at the 
termination shock: the primary interstellar population and the 
secondary population, which comes up between the heliopause 
and the bow shock from charge exchange between the atoms 
from the original interstellar populations and the local heated 
and compressed plasma. These results are presented in Sect. 14.1 
and will be referred to as the composite model, in the sense that 
it is composed of the two populations: primary and secondary. 
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Fig. 3. Classical/Doppler density ratios in the composite model for the solar minimum (left panel) and solar maximum conditions 
(right panel). The Sun is at the center in the (0,0) point, the upwind - downwind direction is along the y = line. The upwind, 
crosswind, and downwind axes are marked with heavy lines. 
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Fig. 4. Density excess q of neutral interstellar hydrogen in the inner heliosphere, defined as a ratio of results of the classical hot 
model to the results of the Doppler model. Left panel: the density ratio at 1 AU as a function of the offset angle from the upwind 
direction, right panel: the density ratio at crosswind as a function of the heliocentric distance. Solid lines correspond to the solar 
minimum conditions, and broken lines to the solar maximum conditions. 



In the second group of simulations, we checked the robustness 
of the conclusions drawn in Sect. 14. II against uncertainties in the 
bulk velocity and temperature of interstellar gas at the termi- 
nation shock and against uncertainties in the ionization rate in 
the inner heliosphere. The results of this series are discussed in 
Section l4~2l 

The calculations were performed on a mesh of heliocentric 
distances r equal to 0.4, 0.7, 1, 2, 4, 7, 10, 20, 50, and 100 AU 
and of offset angles from the upwind directions 9 equal to 0°, 
30°, 60°, 90°, 120°, 150°, and 180°, separately for solar mini- 
mum and solar maximum conditions, characterized by the fol- 
lowing values of the line-integrated solar Lyman-ff flux and net 
ionization rate: 

solar max.: 7 tot = 5.50 x 10 11 photons crrT 2 s _1 ,/3 = 8 x 10 4 s 4 
solar min.: 7 tot = 3.53 x 10" photons crrT 2 s 4 ,/? = 5 x 10 4 s 4 . 

The boundary conditions of the composite simulation were 
the parameters returned by the Moscow MC model at the nose 
of the termination shock for the following set of LIC conditions 
dlzmodenov et alj 12003): proton density « Pi lic = 0.06 cm 4 , 
neutral H density rtH,uc = 0.18 cm 4 , temperature of the gas 
?lic = 6400 K, bulk velocity v^lic = 26.4 km/s, density of He + 
in the LIC n H e+ = 0.008 cm 4 . 



The resulting parameters of the primary and secondary pop- 
ulations at the nose of the termination shock were the following: 

Primary: «TS P n = 0.03465 cm 4 , Vxspri = 28.512 km/s, 
7- T s,pri = 6020K; 

Secondary: «TSsec = 0.06021 cm 4 , vjssec = 18.744 km/s, 
T'TS.sec = 16300 K. 

These parameters were adopted as the boundary conditions 
for the composite simulations, performed with the use of the 
Warsaw model. 

In addition to the dependence of radiation pressure on v r , the 
code also included a simplified model o f the electron ionization 
rate, as adopted by Bzowski et al.l d2008l) . The model is spheri- 
cally symmetric and identical for the solar minimum and maxi- 
mum conditions. Inclusion of this model is a departure from the 
approach adopted in this study to keep all aspects of the model- 
ing as close to the classical hot model as possible. We decided 
to include the electron ionization to make certain that we do not 
miss possible disturbances of the Doppler-to-hot model ratios 
due to the extra ionization term that does not conform with the 
1 jr 1 distance dependence of the charge exchange and photoion- 
ization rates. 

The simulations were repeated with the use of a model that 
differs from the former one only by the lack of sensitivity of 
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Fig. 5. Density ratio q of the classical hot model to the Doppler 
model at the downwind axis as a function of the heliocentric 
distance. Solid line: solar minimum conditions, broken line - 
solar maximum conditions. 



the radiation pressure to vy, i.e. effectively using the classical 
hot model with the electron ionization term added. The radiation 
pressure was equal to the (appropriately scaled) spectral flux av- 
eraged by ±30 km/s about 0, as indicated in Fig[T] The numerical 
accuracy of the solutions was at the level of 1 to 2%. 

The simulations discussed in Sect. 14.21 assume there is only 
one population of interstellar hydrogen at the termination shock. 
We explored the parameter space adopting the net ionization 
rates at 1 AU (3, bulk velocities v sre , and temperatures r src cov- 
ering the entire range of expected values from the low values of 
f3 relevant for the fast solar wind conditions and solar minimum 
to the high values relevant for slow solar wind conditions and 
solar maximum, the ranges of v src and r src from the extremes 
predicted for the primary and secondary populations at the ter- 
mination shock. 

This scan over /3 was performed by assuming that the gas 
temperature and velocity at the termination shock a re equal, re- 
specti vely, to 12500 K and 22 km/s, as derived by ICosta et alj 
d 19991) from an analysis of the heliospheric Lyman-a glow in the 
upwind hemisphere, for 7 tot relevant for solar minimum and max- 
imum conditions(see above). The scan over v src and r src space 
assumed that /3 and I lot values are as listed for solar minimum 
and maximum conditions. 

The baseline set of simulations used the newly-developed 
Doppler model. Comparison simulations used the classical hot 
model on identical spatial mesh and assumed identical param- 
eters at the termination shock and an identical ionization rate. 
In this section the electron ionization was not included, so the 
ionization rate was fell off as 1/r 2 . 

4. Results 

4. 1 . Composite model 

The simulations show that neglecting of the vv-dependence of 
radiation pressure produces an excess of local hydrogen density 
(see Fig. [3). The excess q is defined here as the ratio of density 
returned by the classical hot model n c i aS s to the density returned 
by the Doppler model «dopp : <l = "ciass/ n dopp- Discrepancies be- 
tween the classical and Doppler models appear at ~ 3 AU and 
increase towards the Sun. Depending on the phase of the solar 
cycle and the offset angle from upwind, at 1 AU they range from 
a few dozen percent up to a factor of ~ 1.7 (see Fig. |4j». The ex- 
cess density is a strong function of solar activity; while at solar 
minimum its typical values at 1 AU are about 25%, they are a 



little more than twice as high at solar maximum. An exception is 
a conical region around the downwind axis with the opening an- 
gle of ~ 60°, as shown in Figs.|4]and|5] On the downwind axis, 
the hot model predicts a density deficit, which almost does not 
depend on the activity phase except at the closest distances to the 
Sun. During solar maximum the excess is predicted even in the 
downwind region, though only for distances under ~ 1.5 AU. 

The existence of the density excess can be understood eas- 
ily after inspecting the shape of the solar Lyman-or line profile 
shown in Fig. Q] Since the profile has a minimum close to 
radial velocity, the H atoms before approaching the Sun sense 
a stronger repulsion than in the case when one adopts a "flat" 
radiation pressure with the mean value averaged between +30 
km/s around the line center. As a result, in the Doppler case a 
larger portion of these atoms will be slowed down and repelled 
from the Sun, which reduces the density in comparison with the 
classical case. The ensemble of surviving H atoms has a lower 
bulk velocity than predicted by the classical model, which leads 
to a further enhancement of the excess, because the ionization 
has more time to eliminate the atoms traveling more slowly. 
On the other hand, a sub-population of atoms exists that on ap- 
proach towards the Sun experience a lower radiation pressure 
in the Doppler model than in the classical model. These are the 
atoms from the slow wing of the distribution function. Although 
the Doppler model predicts for them less deceleration than the 
classical model, these atoms contribute comparatively little to 
the entire ensemble because they are more readily eliminated by 
ionization than the faster atoms. 

The Doppler model predicts different ratios of densities be- 
tween pairs of offset angles 9\, 62 (e.g. upwind/crosswind, up- 
wind/downwind etc.) than the predictions of the classical hot 
model. The discrepancies are not strong in the upper hemi- 
sphere, but escalate with the increase in the difference 62 — 61. 
This must be one of the reasons for which interpreting down- 
wind observations in the heliospheric Lyman-a glow has al- 
ways been more challenging than in th e upwind ones and usu- 
ally returned different conclusions (e.g. lLallement et al] [l985a; 
Oue merais et alj Il992l) . The challenge can be better appreci- 
ated when one realizes that the wavelength-sensitivity effects 
are convolved with the effects related to variations in the solar 
Lyman- a flux and of the ion iz ation rate, as discussed, e.g., b y 
Rucinski & Bzowski (Il995bl) : iBzows ki & Rucinskij (Il995albl) : 
iBzowski et al.l (120021) . and with the effects of inhomogeneous 
distribution o f the parameters of inter stellar gas along the termi- 
nation shock (Izmo denov et alJfeOOll) . 

The Doppler and classical models predict surprisingly small 
differences between the local bulk velocities of the gas (Fig. [6}. 
During solar minimum they are about 1 km/s, and they increase 
to ~ 2 km/s during solar maximum. Differences in the radial 
component of the bulk velocities at 1 AU are also small (Fig. 
|7). They begin only at the offset angle 6 60° and increase 
towards downwind, and the classical model systematically pre- 
dicts somewhat higher values. Since the differences in v r are 
small, we do not expect big differences in the model spectral 
profiles of the heliospheric glow. 

Even though the differences in the local velocities are small, 
the effect of the wavelength dependence of radiation pressure on 
the magnitude of fluxes of the two heliospheric neutral H pop- 
ulations is not negligible because of the appreciable differences 
in the local densities of the gas. As shown in Fig. [9] the classical 
excess of the flux is 15 to 20% during solar minimum, depend- 
ing on the population, and during solar maximum it may reach 
a factor of 1.4 to 1.9 in the upwind hemisphere and even more 
in the downwind region, but with the absolute magnitude of the 
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Fig. 6. Bulk velocities of the primary and secondary populations of neutral interstellar H at 1 AU, calculated using the Doppler 
and classical models, shown as a function of the offset angle from upwind for the solar minimum (left panel) and solar maximum 
conditions (right panel). Primary population: solid line (classical) and broken line (Doppler); secondary population: dotted line 
(Doppler) and dash-dot (classical). 
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Fig. 7. Radial velocities of the primary and secondary populations of neutral H atoms at 1 AU as a function of the offset angle 
from upwind for the solar minimum (left panel) and maximum conditions (right panel), calculated with the use of the Doppler 
and classical hot models. The upper pairs of lines in both panels correspond to the primary population and the lower pairs to the 
secondary population. The classical hot model results are drawn, correspondingly, with solid and dash-dot lines, and the Doppler 
model results with broken and dotted lines. 
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Fig. 8. Absolute fluxes of neutral interstellar H atoms at 1 AU from the primary and secondary populations for the solar minimum 
(left panel) and solar maximum conditions (right panel), shown as a function of the offset angle from upwind, calculated with the use 
of the Doppler and classical hot models. The classical model results are drawn with solid line (primary populations) and dash-dotted 
line (secondary population); the Doppler model results with dashed (primary) and dotted lines (secondary). 



flux reduced by 2 orders of magnitude from the solar minimum 
level (Fig. [8]). In the upwind hemisphere it is very weakly sensi- 
tive to the offset angle. The wavelength dependence of radiation 
pressure must then be appropriately taken into account when in- 
terpretating direct in situ observations of neutral interstellar hy- 
drogen atoms, such as those planned for IBEX. The local flux 
of the atoms is a product of the local velocity vector and of the 



local density, which is directly proportional to the density at the 
termination shock. Because of this, the density at the termination 
shock derived from an in situ-measured flux at 1 AU neglecting 
the effect of v r dependence of radiation pressure will be at least 
15% off the mark (Fig|9]l, and most probably more, up to a factor 
of 2. 
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Fig. 10. Deflection of the local upstream direction of the primary and secondary populations of interstellar hydrogen at 1 AU as a 
function of the offset angle from the gas inflow direction (upwind) for solar minimum and maximum conditions. At solar minimum 
(left panel), the primary population is less deflected (solid line: classical model, broken line: Doppler model); the secondary is drawn 
with the dotted (Doppler) and dash-dotted line (classical). Same scheme for the solar maximum plot (right panel). 
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Fig. 9. Classical/Doppler ratios of the fluxes of the primary and 
secondary populations of neutral interstellar H atoms for the so- 
lar minimum and maximum conditions as a function of the offset 
angle from upwind. The lower pair of lines corresponds to solar 
minimum, the upper pair to solar maximum. Within the solar 
minimum (lower) pair, the solid line corresponds to the primary 
population, the broken line to the secondary population. Within 
the solar maximum (upper) pair, the dash-dot line is for the pri- 
mary, the dotted line for the secondary population. 

Another factor that can potentially affect interpretation of 
such measurements are differ ences in the l ocal bulk velocity vec- 
tors of the two populations. Mobi uset al.l 12001) proposed to use 
the relative positions of the beams from the primary and sec- 
ondary populations to check whether their flow directions at the 
termination shock are parallel or not. Such measurements will 
hopefully b e realized soon by a for t hcoming NA SA SMEX mis- 
sion IBEX (iMcComas et alJl2004 120051 120061) . In Fig. Q1 we 
show that the deflections of the two populations differ from the 
predictions of the wavelength-independent model. For both pop- 
ulations in the region of the Earth orbit where detection by IBEX 
is the most probable, the differences in deflection angle are only 
a few degrees and of a similar, but not identical magnitude. The 
magnitude of the deflection is a strong function of solar activity, 
with radiation pressure clearly playing the dominant role, but 
higher values are consistently returned by the Doppler model. 

Even though the differences are just a few degrees, they 
cannot be regarded as negligible. lOuemerais et al. (1999|) and 
lLallem ent et al] d2005l) suggested an upwind direction of neu- 
tral interstellar hydrogen that differs by about 5° from the i nflow 
direction of interstellar helium measur ed by IWittel (12004). The 



difference is commonly attributed to a distortion of the helio- 
spheric interface by the extraheliospheric magnetic field; hence 
mostly affected should be the secondary population and the dif- 
ference in the upwind direction of the primary and secondary 
populations inferred from analysis of an IBEX-like in situ mea- 
surements at 1 AU with the use of a v r -insensitive model might 
be erroneously attributed to a deformation of the interface. Our 
results show that the deflections of the secondary population ob- 
tained from the Doppler and classical models are very similar to 
each other, but the deflections of the primary population differ 
by a few degrees. 

Finally, we assess the influence of the v, -sensitive radiation 
pressure on the expected PUI fluxes at 5 AU crosswind, i.e., 
the locati on where the H + PUI flux observed by Ulysses is the 
strongest. Bzo wski et all (120081) point out that the local produc- 
tion rate of pickup ions, which is a product of the local density 
of neutral hydrogen and of the local ionization rate, is a weak 
function of radiation pressure, hence, necessarily, of its details 
including the v r dependence. They provide an estimate of the 
importance of the v, -dependence of radiation pressure for the 
local production rate of PUI, which suggests it is on the level 
of a few percent. Here we report its influence on the total flux 
of PUI, which was compute d fo llowing the simple ap proach by 
IVasvliunas & Siscodd 19761) and lGloeckler et al.l (11993b . Results 
of these simulations suggest that the classical hot model yields 
an excess flux with respect to the Doppler model, which at 
~ 5 AU crosswind is equal to only ~ 7.5% at solar minimum and 
to 8.5% at solar maximum. Thus we conclude that neglecting the 
v r dependence of radiation pressure results in an overestimate of 
the H + PUI flux at 5 AU from the Sun by about 10%, irrespective 
of the solar activity phase. 

4.2. Comparison with the hot model: a scan of the parameter 
space 

The density excess q is a weakly increasing and almost lin- 
ear function of the ionization rate (Fig. [TTJ; i.e., the higher the 
ionization rate taken, the lower the quality of the approxima- 
tion given by the classical hot model obtained. This can be ex- 
plained by preferential eliminati on of slower atoms from the 
local ensemble. As discussed by lLallement et al.l d!985al) and 
iBzowski et al.1 d!997l) . ionization results in a net acceleration of 
the gas by a few km/s because the atoms that have higher spe- 
cific velocities in the ensemble survive against ionization bet- 
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Fig. 11. Classical hot/Doppler model density ratios for the upwind (dots), crosswind (dashed), and downwind axes (solid) at 1 AU 
as a function of the ionization rate for the solar minimum (left panel) and maximum conditions (right panel). 



ter. Such atoms also have higher radial velocities, so they expe- 
rience a higher radiation pressure due to the non-flat shape of 
the solar Lyman-a line profile. This force repels them from the 
Sun stronger than when the solar line profile is flat. As a result, 
we have a reduction in the density; consequently, the classical 
model predicts an excess of the density, whose magnitude in- 
creases with the increase in the ionization rate. 

The sensitivity of the excess to the ionization rate increases 
with the offset angle from upwind and with the magnitude of 
the solar flux 7 tot . At the upwind axis, the slope dq/d/3 for an 7 tot 
value corresponding to the solar minimum conditions is equal 
to half of the value for crosswind and to ~ 0.1 of the value for 
downwind. For an 7 tot value related to the solar maximum condi- 
tions, the slope is steeper and the highest increase occurs at the 
upwind axis. 

In the downwind region, increasing the ionization rate im- 
proves the quality of the approximation provided by the classical 
hot model with respect to the Doppler model. The parameter q 
is still an increasing function of the ionization rate, but since it 
is lower than 1 in the downwind region, its increase means an 
improvement. 

The magnitude of the excess q is a weak function of the 
temperature and velocity of the gas at the termination shock, as 
shown in Fig. [T2] At 1 AU crosswind, for 7 to t related to the so- 
lar minimum conditions, it varies from 1.15 to 1.30, depending 
on the bulk velocity and temperature of the gas at the termina- 
tion shock; for 7 tot related to the solar maximum conditions the 
amplitude is similar, but the value of q higher. 

Generally, the quality of the hot-model approximation of 
the density is consistently lower for higher values of 7 t ot, which 
is modulated by the solar activity level. The sensitivity of ra- 
diation pressure to the radial velocity thus modifies the re- 
lations between the densities during various phases of solar 
activity (the amplitude of modulation), as discussed, e.g ., by 
Rucin ski & Bzowski ( 1995b); lBzowski & Rucinskil(ll995al ) and 
Bzowski et ail (120021) . 

On the other hand, it seems that the effect on the model val- 
ues of the heliospheric Lyman-o' glow should not be greater than 
10%. An exact calculation requires developing a model, where 
the local bulk velocity and temperature of the gas and their ra- 
dial gradients would be taken into account, together with the 
wavelength dependence of the solar flux and its wavelength- 
differential attenuation, which increases with the increase in the 
solar distance. We believe this effect to be relatively weak be- 
cause during solar minimum, when the contribution to the net 
backscatter intensity from the gas near the Sun is greatest, the q 
factor is lowest. During solar maximum, when q is greatest and 



hence the quality of the v r -independent model is lowest, the lo- 
cal density of the gas is very much reduced, so increased is the 
contribution to the net signal of more distant regions of the he- 
liosphere, where the q values are equal to 1 . In consequence, the 
v r dependence of the radiation pressure should appreciably af- 
fect in situ measurements within ~ 3 to 5 AU from the Sun, but 
the photometric observations of the heliospheric glow only to a 
small degree. 

5. Conclusions 

We compared the density, bulk velocity, and flux of neutral inter- 
stellar hydrogen in the inner heliosphere, calculated using either 
the classical hot model or a newly-developed model with the ra- 
diation pressure being a function of specific radial velocities of 
individual atoms. The conclusions are the following. 

1. Differences between the Doppler and classical hot models 
are restricted to ~ 3 to ~ 5 AU from the Sun, depend- 
ing on the angle off the upwind direction. The classical hot 
model overestimates the density of neutral hydrogen gas ev- 
erywhere in the inner heliosphere apart from a ~ 60° cone 
around the downwind direction, where a density deficit is 
predicted. Generally, the density excess/deficit is a strong 
function of the offset angle from the upwind direction and 
of the heliocentric distance (Figs. [3] @] and|5j. 

2. The magnitude of density excess q varies with the solar activ- 
ity level and is higher at solar maximum (Figs.[3]and|4]i, i.e. 
vv-independent models stronger overestimate the local den- 
sity during solar maximum. 

3. The density excess q is a weakly increasing, almost linear 
function of the ionization rate in the inner heliosphere and a 
weak, almost bilinear function of the bulk velocity and tem- 
perature of the gas at the termination shock (Figs. [TT] and 
H3. 

4. The classical model shows a higher upwind-downwind am- 
plitude of radial velocities at 1 AU than the Doppler model; 
the differences are about 10%. The excess of the absolute 
flux of neutral atoms, returned by the classical model, is 
close to 15% at 1 AU during solar minimum, i.e. at the time 
of observations by the forthcoming IBEX mission. The ex- 
cess is much higher during solar maximum, but the absolute 
values of the fluxes are lower by a few orders of magnitude 
than during solar minimum (Figs. [7|8l9l6t . 

5. The classical model returns a different deflection of the lo- 
cal upstream directions of the primary population of neutral 
interstellar H than the Doppler model. The deflections of the 
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Fig. 12. Density excess q of the classical hot model with respect to the Doppler model of neutral interstellar hydrogen as a function 
of the temperature and bulk velocity in the source region for the solar minimum (left panel) and solar maximum conditions (right 
panel). 



secondary population returned by the classical and Doppler 
models are practically identical, meaning that any interpre- 
tation of in situ measurements of these fluxes must be done 
very carefully, with an appropriate account of the solar line 
profile shape to avoid confusing deflections caused by the 
Doppler effect in the radiation pressure with, e.g., results of 
deformating the heliospheric interface by external magnetic 
field (Fig. [ft). 

6. Discrepancies between the PUI flux at Ulysses in aphelion, 
returned by the Doppler and classical models, are about 10% 
and are almost insensitive to the phase of the solar cycle. 
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